Two-mode theory of vortex stability in multicomponent Bose-Einstein condensates 
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We study the stability and dynamics of vortices in two-species condensates using a two mode model. 
The recent experimental results obtained at JILA (M. R. Matthews, et ai, Phys. Rev. Lett. 83 
(1999) 2498) and theoretical predictions based on the Gross-Pitaevskii equations are verified. We also 
make an exhaustive analysis of the stability properties of the system when the relative populations 
of the two species and/or their relative scattering lengths are changed and prove that stabilization 
of otherwise unstable configurations can be attained by controlling the relative population of both 
species. 
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I. INTRODUCTION 

Vortices appear in many different physical contexts rang- 
ing from classical phenomena such as fluid mechanics Jl| and 
nonlinear Optics [HI to purely quantum phenomena such as 
superconductivity |3j and superfluidity yg. 

In fact a vortex is the simplest topological defect one can 
construct (5) : in a closed path around a vortex the phase of 
the involved field undergoes a 2n winding and stabilizes a zero 
value of the field placed in what is called the vortex core. The 
vortex is usually stabilized by topological constraints since re- 
moving the phase singularity implies an effect on the bound- 
aries of the system which is difficult to achieve using only local 
perturbations. 

In particular, the concept of a vortex is central to our un- 
derstanding of superfluidity and quantized flow. This is why 
after the experimental realization of Bose-Einstein conden- 
sates (BEC) with ultracold atomic gases |J the question of 
whether atomic BEC's are superfluids has triggered the anal- 
ysis of vortices. The main goals have been to propose a robust 
mechanism to generate vortices |7| and detect them Fi- 
nally another research area has been the analysis of vortex 
stability 

Although most of the theoretical effort has concentrated on 
single condensate systems, the final experimental realization 
of BEC vortices [^3j following the proposal of Ref. jl4| was 
attained with a two-species Rb condensate. The two species 
correspond to two different hyperfine levels of Rb, denoted |1) 
and 1 2). In Ref. Q it was found that the vortex is stable in 
only one of the two possible configurations. The stable con- 
figuration corresponds to the vortex placed on the |1) state, 
which is the one with the larger scattering length. The other 
possibility, i. e. when the vortex is placed in the |2) state, 
which is the one with the lowest self-interaction coefficient, 
leads to some kind of instability. To simplify the notation in 
what follows we will use a shorthand notation for these states 
calling 1 1,0) to the state where the vortex is put in |1) and 
|0, 1) to the state where the vortex is placed in |2). 

In a recent work [^5) we have shown, using numerical sim- 
ulations, that the origin of the instability of state |0, 1) is 
purely dynamical and can be understood within the frame- 
work of mean field theories for the double condensate system. 
A consequence of that analysis is to clarify the instability 



mechanism which does not lead to expulsion of the vortex 
from condensate, but to the stablishment of a complex state 
where the phase singularity is periodically trasfered from one 
specie to the other. 

In this paper our intention is to understand in more detail 
the instability mechanism and to study configurations which 
can be obtained experimentally by varying the relative pop- 
ulations of the two species. The results are also applicable to 
other multiple condensate systems where the different scatter- 
ing lenghts have values different from those of Rubidium. To 
make the analysis simpler and to ease the physical interpreta- 
tion we make most of the analysis using a simplified two-mode 
model whose predictions are in close agreement with the ex- 
periment of Ref. [^) and the numerics of Ref. [JH| as well as 
with numerical simulations of the full three-dimensional mean 
field equations ruling the phenomenon. 

Our plan is as follows: In Sec. [ij] we present our problem in 
appropriate form and obtain the simplified equations for the 
two-mode model. In Sec. [ill] we apply the two-mode model 
to the situation discussed in the previous analysis In 
Sec. [iy| we make a complete stability analysis of the rele- 
vant configurations and make some predictions which are ex- 
perimentally testable. Finally in Sec. |v] we summarize our 
conclusions. 



II. THE MODEL 

A. Mean field equations for the two-condensate 
system 

In this work we will use the zero temperature approxi- 
mation, in which collisions between the condensed and non 
condensed atomic clouds are neglected. In the two species 
case this leads to a pair of coupled Gross-Pitaevskii equations 
(GPE) for the condensate wavefunctions of each specie 
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Where Uij = 4irh 2 aij/m, a%j being the s-wave scattering 
lenghts of 1 — 1, (an), 2 — 2 (022), and 1 — 2 (012) binary 
collisions. 

To simplify the formalism we assume that both potentials 
are of the form V\{r) — V2(r) = ^mco 2 r 2 . 

Next we change to a new set of units based on the trap 
characteristic length, ao = yj%/mu, and period, r = 1/w 
defined as x — > x/ao, t — > t/r, Uij = 4nNja,ij/ao and ^(x) = 
Njipj(x). Equations ([j]) conserve the number of particles on 
each hiperfine level so that we may choose 



= / \M^)\ = 1. 
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This choice implies that the particle number of each specie 
appears on the nonlinear coefficients Uij. 

The experimental results Jl^] and our previous theoretical 
analysis |l5| correspond to systems in which the number of 
particles is the same for each component, Ni = N2 — N, 
but in general one could allow any proportion between the 
populations of the different levels. 

With the previous reescaling, the GPE for the multicom- 
ponent system read 
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Since the realistic values of Rb scattering lengths are in 
the proportion an : a 12 : a 22 = 1.00 : 0.97 : 0.94 @, the 
coefficients of the matrix of nonlinear coefficients satisfy the 
relations Uii/«i2 = an iVi/ai2-/V2, U21/U22 = a,i2Ni/a,22N2, 
which means that except for the particular case in which 
N\ = N2 = N, this matrix is nonsymmetric. In terms of 
the population imbalace /3 = N2 /Ni , and for a fixed total N 
the matrix can be written as 
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B. Derivation of a two-mode model 

In our previous work jl5| we considered the full problem 
and worked on the basis of the full GPE to prove the insta- 
bility of the stationary solution 0, 1), as well as the stability 
of state 1 1,0) for typical experimental conditions. The sta- 
bility analysis based on the full GPE demonstrated that the 
instability was mediated by the growth of a core mode. This 
fact makes plausible the description of the two-condensate 
dynamics by the use of only two modes for each level one cor- 
responding to the vortex and other to the core mode. This 
approach, which corresponds to retaining the stationary plus 
active modes and has been used succesully in the analysis 
of other nonlinear problems |n| , should work at least in the 
linear regime where the perturbations are small. Mathemati- 
cally, the idea is to approximate 



V>i(x) ~ a(t)^ 9 i(x) + 6(t)V' e i(x), 
V> 2 (x) ~ c(t)V> s2 (x) + d(t)i> e2 (x). 
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where ip g j(x) is the spatial wavefunction of the ground state 
or core mode for species \j) and tpej {%) corresponds to a repre- 
sentation of the single vortex wavefunction. This approxima- 
tion implies some loss of information on the dynamics but is 
not essential for our results as will be shown later. We will see 
that this ansatz reflects the essentials of the dynamics with 
good accuracy. 

tp g , ip e can be choosen as any approximation to the ground 
and first excited states of the single specie equations, provided 
they are orthogonal. Our choice will be to use the eigenfunc- 
tions of the d— dimensional harmonic oscillator which are the 
exact solutions in the linear case and allow simple manipula- 
tion since their analytic form is known: 
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Other choices are possible with the only change of several 
coefficients related to integrals involving ip g and tp e . In our 
treatment we will consider simultaneously the two and three 
dimensional configurations. It was shown in a previous work 
concerning single condensates jllj, that the transition from a 
spherical trap to a pancake preserves the shape and number 
of unstable modes. This fact motivated us in Ref. to use a 
simpler model in which the dependency on the axial direction 
has been integrated out. Our present analysis applies equally 
to the simplified 2D situation used in Ref. jl5| as well as to 
the full 3D problem and proves that there are no essential 
differences between the two and three-dimensional models for 
the phenomena studied here. 

Inserting ansatz (^) into GP equation (^) and proyecting 
on ipgj and ip e j the following set of coupled nonlinear ODE 
are obtained 
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where 71 = (^,^9), 72 = (^ 2 ,^ 2 ),73 = {i>i,i>l) and 
Eg = (ij} g ,Hoi>g),E e = {4> e ,Ho4> e ) being H = -5A + \r 2 . 
For our particular choice of ip g and ip e the numerical val- 
ues of these coefficients are 71 = 1/27T, 72 = 1/47T, 73 = 
l/47r,£; 2D = l,E 2D = 2, 7l M = l/(2rr) 3 / 2 ,72 3ti = 1/2(2tt) 3 / 2 , 
5/12(27r) 3/2 , E z g D = 3/2,E't D = 5/2. 
Equations (pi) satisfy discrete conservation laws correspond- 
ing to the number of particles of each species and angular 
momentum 



73 



14" + |&| =1, 
|c| 2 + |d| 2 = l, 
|6| 2 + |d| 2 = L , 



(8a) 
(8b) 
(8c) 



2 



I/O being the angular momentum of the initial data. There is 
another conservation law for energy which is not relevant for 
our purpouses. 

It is convenient to change to the modulus-phase represen- 
tation as, a = p a e l ^ a ,b — pbe l ^ b , 
([?]) become 



p c e^,d = p d e^. Eqs. 
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Despite the apparent complexity of this system it can be seen 
that the four phase variables can be reduced to only one defin- 
ing $ = 4> b + 0c — 0a — 4>d, and then 
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We have now five equations in ([To]) plus four conservation 
laws, which means that the system can be (at least formally) 
integrated. This fact excludes the possibility of chaotic be- 
havior in the system. 

These equations can be further simplified by defining den- 
sity variables related to pj and X = papbpcpd and using the 
conservation laws. We will not follow this route since all the 
simplified models such as the one presented in Eqs. ( Jlo| ) have 
singularities when any of the densities is zero. This fact makes 
the equations unuseful for the purpouse of stability analysis, 
since the stationary states are singular solutions of these sys- 
tems. 



in the case where the populations of both hiperfine levels 
are equal. It means that /3 = 1 and thus it2i = tti2. It 
is interesting to study the dynamics of small perturbations 
of the initial data a(0) = 0,6(0) = l,c(0) = l,d(0) = 0, 
which physically corresponds to the stationary state |1, 0) and 
a(0) = 1,6(0) = 0,c(0) = 0,d(0) = 1 which corresponds to 
1 1 , 0) . These initial data correspond to periodic solutions of 
the amplitude equations (m), which are 
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respectively. To have a clear picture of what is going on 
we have first simulated the dynamics of these states when 
small perturbations are added, the results being summarized 
in Figs. to ^ 
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FIG. 1. Stability of the configuration 1 1 , 0) . Snapshots of 
the spatial density of (a) |1) and (b) |2). Evolution the am- 
plitudes of the modes with time (c) \a\ (solid line) and |6| 
(dashed line); (d) \c\ (solid line) and \d\ (dashed line). 



III. DYNAMICS IN THE PHYSICALLY 
RELEVANT CASE 

The configurations studied experimentally in ^| and nu- 
merically in |l5j correspond to single vortices in |1) or |2) 
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FIG. 2. Snapshots of the evolution of an unstable vortex 
(state 1 0, 1 } ) - Evolution of the spatial density of (a) 1 1) and 
(b)|2>. 



system while a phase singularity appears in 1 1} and occupies 
the center of the atomic cloud [Ell . This dynamics is recurrent 
as can be seen from the evolution of the relevant variables 
(Fig. 




3 2 10 12 3 
FIG. 3. Evolution of the position of the phase singularity 
corresponding to the simulation shown in Fig. ^. (a) Phase 
singularity in |1) (b) Phase singularity in |2). 
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It is clear in Fig. [I] that the configuration which has a 
vortex in |1) is dynamically stable. However, when the vortex 
is placed in |2) an instability develops and the response to 
small perturbations is to trasfer the vortex to |1) and start 
a periodic trasfer dynamics. The snapshots of the trasfer 

(evolution of the density) are shown in Fig. ^ In Fig. || bers uniformly distributed between and 0.02. rj are random 
it is seen how the phase singularity in |2) spirals out of the 



FIG. 4. Evolution of the amplitudes of the modes (a) 
a | (dashed line) and |6| (solid line); (b) |c| (dashed 
line) and |d| (solid line). The simulation is done with 
a random perturbation to the co nfigura tion with vor- 
i. e., a(0) 

ir4 . 



tex in 1 2) , 
c(0) 
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; d(0) = v/T 



ei and 62 are random num- 



numbers uniformly distributed between and 1. 
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IV. STABILITY THEORY 

A. Problem statement 

Our numerical simulations of the reduced system (Q) show 
that in the equal population case, Ni — N2 — N, and for arbi- 
trary nonlinearities one of the possible stationary states of the 
system is stable while the other is unstable. It is our purpouse 
in this section to make a complete analysis of the stability of 
the system for any proportion of the populations j3 = N1/N2 
and any value of the nonlinear coefficients (e.g. total num- 
ber of particles, N, and scattering lengths Oij). These results 
could be specially relevant to predict for a specific multiple- 
condensate system the existence of stable vortex states. For 
the case of Rb the results can be applied to study the possi- 
bility of stabilizing different configurations. 



B. Stability of state |1,0) 

When a vortex is placed at 1 1} the resul ting stationary state 
is a periodic orbit described by Eq. (11a). Its direct stability 
analysis using Eqs. (j7) would lead to time-dependent pertur- 
bation equations which should be analyzed using Floquet the- 
ory. A way to circunvent partially this situation is to change 
to the rotating system defined as 
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and the stationary solution is an equilibrium point of Eqs. 
(fili|): 5o = 0,&o = l,co = l,do = 0. To study its stability 
we linearize Eqs ( |l3| ) around the equilibrium point and define 
the perturbations through 



5 = 5o + 5 a (t), 
b = b +6 b (t), 
c = c + S c (t), 
d = d + 6 d {t). 
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where A a = 112271 + 112172 — 272W11 — 711*12, A_ d = U1173 + 
111272 — 2ii2272 — 112173- The perturbations for 6 and c have 
a neutral behavior and their evolution is ruled by quadratic 
terms. If we write the equations for the perturbations and 
their complex conjugates to obtain the full stability spectrum 
we have 
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The eigenvalues of this matrix can be obtained analytically, 
the result being 
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There is only one stability condition which is |A a + A d | > 
272v^2iWi2- Since the Aj are functions of 7,-, which depend 
on the dimensionality, we must distinguish now results for the 
2D and 3D cases separately. In the two dimensional case we 
obtain 



Mil + U12 > 2v / U2lWl2, 

while for 3D the condition is 
7 1 

771*11 — 77W21 + «12 > 2^/ 11211112- 
6 D 
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Taking into account the fact that the numerical values of an 
and ai2 are very close we find that Eqs. (|l^) and ([il]) are very 
simmilar. This is why we will use one of them, Eq. (|l8|), for 
the subsequent analysis. If we write these equations in terms 
of P and use the scattering length values of Rb, we obtain the 
following stability condition 



oil 
ai2 



+ /3>2v^. 
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For Rb values the inequality (|20|) is allways satisfied, which 
proves that the configuration with a vortex in |1) is always 
linearly stable no matter what is the relation between the pop- 
ulations P, which implies also the stability of the experimental 
configuration with vortex in |1) studied in Refs. |l3[|lE| , 

It is also interesting that the stability properties do not de- 
pend on the total number of particles but only on the relation 
between the populations |22|. 



Their evolution laws are 
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C. Stability of the vortex in 2) state 

The stability analysis of configuration |0, 1), which corre- 
sponds to initial data ao = 1, 60 = 0, Co = 0, do = 1, is com- 
pletely equivalent to the previous one. In fact, elementary 
symmetry arguments allow to ensure that the result should 
be equivalent with the only interchange of the indices 1^2, 
i.e., the stability condition is 

^lp + l>2^ (21) 

(221 

This inequality is not verified out for a range of (5 values. 
Solving the algebraic equation for f3 one finds for the critical 
values 




For the case of Rb it is found that the unstable range is 
P £ [0.73, 1.49]. This result is interesting since it means that 
there are choices of the population imbalance /3 which allow 
stabilization of the vortex in |2). We have analyzed the ra- 
tio of angular momentum trasfer from component |2) to 11} 
as a function of f3 from numerical simulations of Eqs. (pj). 
The results are presented in Fig. [|. This is one of the main 
results of the paper and a prediction which can be experi- 
mentally tested. We have also made compared these results 
with numerical simulations of the Gross-Pitaevskii equation 
with good qualitative agreement. Of course the exact range 
depends on the size of the perturbation added to destabilize 
the vortex, this means that Fig. 5 must be taken only as ori- 
cntative since the exact profile can change when the size of 
the perturbations added are changed. However the f3 stability 
ranges found by numerical simulations for different perturba- 
tions are simmilar to the prediction of the two mode theory. 




0.2 0.6 1.0 1.4 1.8 2.2 

P 

FIG. 5. Angular momentum trasfer as a function of f3. All 
the angular momentum is put initially at component |2). (a) 
We simulate the dynamics and plot the maximum fraction of 
angular momentum at component |1), /, which is a measure of 
the instability, as a function of /3. The vertical lines mark the 
points where the stability analysis predicts instability of the 
configurations, (b) Value of the real part of the eigenvalues 
leading to instability A = iirNanq/ao (c) Results of realistic 
numerical simulations of the fraction of angular momentum 
trasfered as a function of . 



D. Nonlinear stability analysis 

The nonlinear stability analysis is quite technical and dif- 
ficult for this particular case and will be the subject of future 
work. However, there are indications that the stable configu- 
ration 1 1, 0} should be sensitive to appropriate small finite am- 
plitude perturbations. The reason is very simple, if one looks 
to Fig. 4(b) when the vortex is almost completely trasfered 
to state 1 1) (for a time near 20 time units) we are in a config- 
uration which is close to the stable state |1,0), but which is 
unstable. 

In fact we have added finite amplitude perturbations to the 
configuration |1,0) and find that a periodic trafer dynamics 
is also induced very simmilar to the dynamics of |2, 0). The 
main difference is that |1,0) is linearly stable, which makes 
this configuration more robust but yet not completely stable. 

V. CONCLUSIONS AND DISCUSSION 

Summarizing the work presented in this paper we have 
completed the task of analyzing stability properties of vortices 
in double Rb condensates, although our theory is much more 
general and can be applied to any two (multiple)-condensate 
system. Our treatment is based on a simplified two-mode 
model which captures many of the relevant dynamical fea- 
tures of the problem. We have attained several goals in this 
work. First, the stability results of Refs. Jl^,[l5| are repro- 
duced for the Ni = N2 case. Second, the instability mecha- 
nism consisting in vortex exchange between the two species is 
supported and described in detail here. Third, we rise a new 
prediction which consists on the fact that population imbal- 
ances can stabilize vortices in |2) states and also prove that 
vortices in |1) can be destabilized by adding finite amplitude 
perturbations to the initial data. These predictions can be 
tested with current experimental setups and can be another 
tests of the existence of purely dynamical instabilites in the 
two-condensate system of Ref. |Tl3| . 

The proposed possibility of making the condensate in |2) 
stable or unstable by controlling the population ratio is in- 
teresting from the viewpoint of condensate engineering. We 
hope that this works helps in the task of understanding the 
complex dynamics of vortices in Bose-Einstein condensates. 
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